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ABSTRACT 

During the past thirty years, numerical methods 
and simulation tools for incompressible flows 
have been advanced as a subset of CFD discipline. 
Even though incompressible flows are 
encountered in many areas of engineering, 
simulation of compressible flow has been the 
major driver for developing computational 
algorithms and tools. This is probably due to 
rather stringent requirements for predicting 
aerodynamic performance characteristics of flight 
vehicles, while flow devices involving low speed 
or incompressible flow could be reasonably well 
designed without resorting to accurate numerical 
simulations. As flow devices are required to be 
more sophisticated and highly efficient, CFD tools 
become indispensable in fluid engineering for 
incompressible and low speed flow. This paper is 
intended to review some of the successes made 
possible by advances in computational 
technologies during the same period, and discuss 
some of the current challenges. 

I. INTRODUCTION 

Almost all, if not all, fluid is compressible. 
Incompressible flow is an approximation of flow 
where flow speed is insignificant everywhere 
compared to the speed of sound of the medium. 
Mathematical expressions of incompressible flow, 
such as the incompressible Navier-Stokes 
equations, are therefore approximations to the real 
world. If incompressible flow is defined this way, 
the majority of the fluid and associated flow we 
encounter in our daily lives belong to the 
incompressible category. For example, most of 
the flow associated with air and water, (such as 
flow related to automobile, ships, submarines, 
water supply through pipes and channels, 
hydraulic turbine, pumps, low speed airplane and 


trout swimming in mountain stream) and flow of 
biofluid (such as blood) are all in the 
incompressible flow domain. One of the earliest 
mathematical models of incompressible flow is 
the famous equation by Bernoulli, who developed 
the model equation while investigating blood flow 
in 1730 (see an interesting article by Quarteroni 
[1]). It is not surprising to find many attempts to 
investigate incompressible flow both analytically, 
experimentally and computationally ever since. 

Computational study of incompressible flow 
problems both in basic research and in industrial 
applications has been performed for several 
- decades. Numerical solutions for such basic fluid 
dynamics problems as flow past a circular 
cylinder, flow through channel, duct and pipe, 
flow over a backward facing step were presented 
as early as 1930s (e.g. Thom [2] for circular 
cylinder). In the computational fluid dynamics 
(CFD) community, especially in aerospace, CFD 
has been synonymous to computational 
aerodynamics. Computational incompressible 
fluid dynamics was treated as a special case of 
aerodynamics, that is, the low speed limit of 
compressible fluid dynamics. This is probably due 
to rather stringent requirements for predicting 
aerodynamic performance characteristics of flight 
vehicles, while flow devices involving low speed 
flow can be reasonably well designed without 
resorting to accurate numerical simulations (e.g. 
hydraulic turbine for hydroelectric power plants 
and liquid fuel and oxidizer devices for the Space 
Shuttle main engine were designed without help 
from CFD approach). As flow devices become 
increasingly compact and highly efficient, 
pushing the conventional operating envelope, 
requirements on incompressible CFD tools 
became more demanding just as aerodynamic 
performance prediction tools require quantitative 
prediction capability. This trend is reflected into 
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the development of various incompressible flow 
solution methods and associated tools, especially, 
in conjunction with high-fidelity computations 
using high-end computing facilities. 

Mathematically, the incompressible flow 
formulation poses unique issues not present in 
compressible equations because of the 
incompressibility requirement. Physically, 
information travels at infinite speed in an 
incompressible medium, which imposes stringent 
requirement to computational algorithm of 
satisfying incompressibility as well as difficulties 
in designing downstream boundary conditions. 
Differences in various methods of solving the 
incompressible flow equations originate from 
different strategies of satisfying incompressibility. 

There are generally two different approaches in 
solving the incompressible flow equations, i.e. 
incompressible Navier-Stokes equations. The first 
one is based on satisfying the incompressibility 
directly. If primitive variables, pressure and 
velocity, are chosen, pressure is used as a 
mapping parameter to satisfy incompressibility. 
This class of methods is generally known as 
“pressure-based” method. Instead of pressure and 
velocity, derived quantities like stream fiinction- 
vorticity and vorticity-velocity can be used 
resulting in different sets of governing equations. 
For general three-dimensional applications, 
however, the primitive variable formulation poses 
the least difficulty in geometry modeling and in 
setting the boundary conditions. The second 
approach is based on compressible flow 
formulation where momentum and continuity 
equations are coupled through the use of density. 
Incompressibility is recovered as a limiting case 
of this formulation. This class of methods is 
known as “density-based” method. The artificial 
compressibility method is a typical of this 
approach, which is a special case of 
preconditioned compressible flow methods. 

During the last several years, a large number of 
review articles and books on CFD discussed 
incompressible flow methods. For more 
comprehensive reviews of computational methods 
for incompressible flow in general, readers are 
referred to these materials, i.e. Roach [3], Peyret 
and Taylor [4], Hirsch [5], Kwak [6], Gunzburger 
and Nicolades [7], Hafez and Oshima [8], Gresho 
and Sani [9] and Hafez [10]. In these books and 
articles, formulations, numerical methods and 
solutions to fundamental fluid dynamic problems 
are presented fairly extensively. In the present 


paper, a short summary of methods from a 
historical perspective is given. After thirty some 
years of numerical simulation of incompressible 
flows, industrial problems involving complex 
systems are now solved rather routinely. The 
current review is focused on these applications to 
answer questions like what have been 
accomplished to date and what challenges still 
remain to be resolved. There are a vast number of 
these cases, and those presented here are, of 
course, only limited samples to illustrate the 
author’s opinion. 

II. EVOLUTION OF SOLUTION METHODS 

For general applications, the formulation based on 
primitive variables is the most flexible choice. 
Other formulations using derived variables, such 
as stream function and vorticity, eliminate the 
need for pressure boundary conditions, however, 
at the expense of adding a new requirement for 
derived variables such as vorticity boundary 
conditions. All these formulations produce 
comparable results for two-dimensional problems. 
However, extension to three dimensions is not 
straightforward for stream function formulation. 
For vorticity-velocity formulation in three 
dimensions, three Poisson-type vorticity equations 
are introduced that makes the formulation very 
expensive to solve. In the present paper, examples 
of three-dimensional applications are primarily 
obtained using primitive variable formulation. 
Therefore, in this section, evolution of the two 
most popular approaches using primitive variables 
is briefly reviewed. Other formulations and 
examples can be found in the literature [3-10], 

2.1 Formulation 

Three-dimensional incompressible flow with 
constant density is governed by the following 
Navier-Stokes equations: 


ofy du i u ] dp dx^ 
dt dxj dx i dxj ( ' 2 ' 2 ^ 

Where t is the time, x* the Cartesian coordinates, 
Ui, the corresponding velocity components, p the 
pressure, and the viscous stress tensor. Here, 
all variables are non-dimensionalized by a 
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reference velocity and length scale. For turbulent 
flow, T- should include turbulent stresses. 

The major difference between the incompressible 
and the compressible Navier-Stokes formulation 
is in the continuity equation. The incompressible 
formulation can be viewed as a singular limit of 
the compressible one. Satisfying the mass 
conservation equation, therefore, is the primary 
issue in solving the above set of equations. 
Physically, incompressible flow is characterized 
by elliptic behavior of the pressure waves, the 
speed of which in a truly incompressible flow is 
infinite. The pressure field is desired as a part of 
the solution. However, the pressure condition is 
needed on the boundary for numerical 
computation. One can utilize other formulation 
using derived quantities like vorticity to eliminate 
the pressure from the boundary conditions; 
however, boundary conditions for the derived 
variables have to be imposed. In realistic 3-D 
problems, these derived quantities are difficult to 
define or impractical to use. The primitive 
variable formulation, namely, using pressure and 
velocities as dependent variables, then becomes 
more convenient and flexible in 3-D applications. 
However, in this formulation, mass conservation 
and its relation to pressure has to be properly 
handled while achieving computational efficiency. 
Various techniques have been developed in the 
past, none of which have been proven to be 
universally better than the others. 

2.2 Methods Based on Pressure Iteration 


In 1965, Harlow and Welch [11] published the 
first primitive variable method using a derived 
Poisson equation for pressure. In this method, 
called the marker-and-cell (MAC) method, the 
pressure is used as a mapping parameter to satisfy 
the continuity equation. By taking the divergence 
of the momentum equation, the Poisson equation 
for pressure is obtained: 


2 3hi d du. 
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Where 


du i u j dr tj 
dx , dx { 


The usual computational procedure involves 
choosing the pressure field at the current time step 
such that continuity is satisfied at the next time 


step. The original MAC method is based on a 
staggered arrangement on a 2-D Cartesian grid. 
The staggered grid conserves mass, momentum, 
and kinetic energy in a natural way and avoids 
odd-even point decoupling of the pressure 
encountered in a regular grid [12]. However, these 
properties become unclear when applied to 
generalized curvilinear coordinates. Even though 
the original method used an explicit Euler solver, 
various time advancing schemes can be 
implemented here. Ever since its introduction, 
numerous variations of the MAC method have 
been devised and successful computations have 
been made. 

The major drawback of this method is the large 
amount of computing time required for solving 
the Poisson equation for pressure. When the 
physical problem requires a very small time step, 
the penalty paid for an iterative solution procedure 
for the pressure may be tolerable. But the method 
as a whole is slow and the pressure boundary 
condition is difficult to specify. One important 
aspect of the numerical solution of the Poisson 
equation for pressure is tied to the spatial 
differencing of the second derivatives. To satisfy 
the mass conservation in grid space, the difference 
form of the second derivative in the Poisson 
equation has to be constructed consistently with 
the discretized momentum equation (see Kwak [6] 
for more detail). 

Since, for a steady-state solution, the correct 
pressure field is desired only when the solution is 
converged, the iteration procedure for the pressure 
can be simplified such that it requires only a few 
iteration at each time step. The best-known 
method using this approach is the Semi-Implicit 
Method for Pressure-Linked Equations (SIMPLE) 
[13-14]. The unique feature of this method is the 
simple way of estimating the velocity correction 
and thus the pressure correction. This feature 
simplifies the computation but introduces 
empiricism into the method. Despite its 
empiricism, the method has been used 
successfully for many computations. This method 
has been the basis for many pressure-based 
methods developed ever since (e.g. [15] and many 
others). It is not, however, the intention of this 
paper to fully review the evolution of this method. 
Many successful applications and many tools 
developed to date are probably the real testimonial 
of the value of this pressure-based approach in 
general. 
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Fractional-step procedure for pressure-based 
method 

The MAC method can be viewed as a special case 
of the projection method (i.e. Chorin, [16]) or, 
from the operator point of a view, one variation of 
the fractional step method [17-19]). Here, 
operator splitting can be accomplished in several 
ways by treating the momentum equation as a 
combination of convection, pressure, and viscous 
terms. The common application of this method is 
done in two steps. The first step is to solve for an 
auxiliary velocity field using the momentum 
equation in which the pressure-gradient term can 
be computed from the pressure in the previous 
time step or can be excluded entirely. In the 
second step, the pressure is computed which can 
map the auxiliary velocity onto a divergence-free 
velocity field. This procedure was first applied in 
staggered Cartesian grid by Harlow and Welch 
[11], and later extended to regular grids (e.g. 
Dwyer [20] and many others). 

One particular aspect of the fractional step method 
requiring special care is the intermediate 
boundary condition. Orszag et al. [21] discussed 
this extensively. Rosenfeld et al. [22] devised a 
generalized scheme where physical boundary 
conditions can be used at intermediate steps. As 
with other pressure-based methods, the efficiency 
of the fractional step method depends on the 
Poisson solver. A multi-grid acceleration, which 
is physically consistent with the elliptic field, is 
one possible avenue to enhance the computational 
efficiency. In light of successful computations in 
Cartesian coordinates using staggered grid, it 
became of interest to develop a method for 
general three-dimensional coordinates. However, 
extension of the staggered arrangement to general 
computational cells is not straightforward and 
involves features related to geometric 
conservation. A general numerical method based 
on a staggered grid was developed at NASA 
Ames among others, initially by Rosenfeld et al. 
[22]. Later, Kiris and Kwak [23] developed a 
more robust implicit procedure for not so “nice 
grid” using finite-volume framework on a 
staggered grid. A staggered grid has favorable 
properties in Cartesian coordinates, such as 
coupling odd-even points. However, it is 
debatable whether a staggered grid has clear 
advantages over a regular grid in generalized 
curvilinear coordinates. In practice, this is 
because that flux balancing in general curved cells 
requires interpolation that makes the differences 
between cell-centered and staggered arrangement 


insignificant. Further discussion related to this 
approach can be found in recent collection of 
papers by Hafez [10] (for example, in papers by 
Loner et. al [24] and more). 


2.3 Artificial Compressibility Method 


Advances in CFD in aerospace have been made in 
conjunction with compressible flow computations. 
Therefore, it is of significant interest to be able to 
use some of these compressible flow algorithms. 
One of the early attempts to do this is the artificial 
compressibility method by Chorin (1967) [25]. In 
this formulation, the continuity equation is 
modified by adding a time-derivative of the 
pressure term resulting in 


1 dp | du i _ Q 
/? dt dx t 


(2.4) 


where (3 is an artificial compressibility or a 


pseudo-compressibility parameter. The two terms, 
artificial compressibility and pseudo- 
compressibility, are used interchangeably in the 
literature. Together with the unsteady momentum 
equations, this forms a hyperbolic-parabolic type 
of time-dependent system of equations. Thus, the 
implicit schemes developed for compressible 
flows can be implemented. It is to be noted that t 
no longer represents a true physical time in this 
formulation. 


Physically, this means that waves of finite speed 
are introduced into the incompressible flow field 
as a medium to distribute the pressure. For a 
truly incompressible flow, the wave speed is 
infinite, whereas the speed of propagation of these 
pseudo waves depends on the magnitude of the 
artificial compressibility. In a true incompressible 
flow, the pressure field is affected instantaneously 
by any disturbances in the flow field, but with 
artificial compressibility, there will be a time lag 
between the flow disturbance and its effect on the 
pressure field. Ideally, the value of the artificial 
compressibility is chosen to be as high as the 
particular choice of algorithm allows, so that the 
incompressibility is recovered quickly. This has 
to be done without lessening the accuracy and the 
stability property of the numerical method 
implemented. On the other hand, if the artificial 
compressibility is chosen such that these waves 
travel too slowly, then the variation of the 
pressure field accompanying these waves is very 
slow. This will interfere with the timely 
development of the viscous boundary layer. In 
viscous flows, the behavior of the boundary layer 
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is very sensitive to the streamwise pressure 
gradient, especially when the boundary layer is 
separated. If separation is present, a pressure wave 
traveling with finite speed will cause a change in 
the local pressure gradient that will affect the 
location of the flow separation. This change in 
separated flow will feed back to the pressure field, 
possibly preventing convergence to a steady state. 
Especially for internal flow, the viscous effect is 
important for the entire flow field and the 
interaction between the pseudo pressure-waves 
and the viscous flow field becomes very 
important. 

Artificial compressibility relaxes the strict 
requirement of satisfying mass conservation in 
each step. However, to utilize this convenient 
feature, it is essential to understand the nature of 
artificial compressibility both physically and 
mathematically. Chang and Kwak (1984) [26] 
reported details of artificial compressibility, and 
suggested some useful guidelines for choosing the 
artificial compressibility parameter. 

For obtaining steady state solutions, the governing 
equations are iterated in pseudo time until the 
solution converges to steady state, at which time 
the original incompressible Navier-Stokes 
equations are satisfied. Time dependent 
calculation of incompressible flows are especially 
time consuming due to the elliptic nature of the 
governing equations. Physically, this means that 
any local change in the flow has to be felt by the 
entire flow field. Numerically, this means that in 
each time step, the pressure field has to go 
through one complete steady-state iteration cycle, 
for example, by Poisson-solver-type pressure 
iteration or artificial compressibility iteration 
method. In transient flow, the physical time step 
has to be small and consequently the change in the 
flow field may be small. In this situation, the 
number of iterations in each time step for getting a 
divergence-free flow field may not be as high as 
regular steady-state computations. However, 
time-accurate computations are in general 
extremely time-consuming. Therefore, it is 
particularly desirable to develop computationally 
efficient methods by implementing a fast 
algorithm and by utilizing computer 
characteristics such as parallel processing. 

Various applications, which evolved from this 
concept, have been reported for obtaining steady- 
state solutions (e.g. Kwak et al. [27], Chang et al. 
[28], Choi and Merkle [29]). To obtain time- 
dependent solutions using this method, an 


iterative procedure can be applied in each physical 
time step such that the continuity equation is 
satisfied. Merkle and Athavale [30] and Rogers et. 
al [31] reported successful computations using 
this pseudo time iteration approach. 

Pre-conditioning 

In incompressible flow, the density is constant, 
which is the source of difficulty in solving the 
compressible Navier-Stokes equation as the flow 
approaches the incompressible limit. Adding a 
non-physical pressure term made it possible to 
adopt a time-like iteration scheme in solving the 
governing equations. The success of this approach 
prompted a series of attempts to develop 
preconditioning methods to be used in 
conjunction with compressible Navier-Stokes 
equations (Turkel [32], Van Leer et. al [33], Choi 
and Merkle [34], Venkataswaran and Merkle 
[35]). This method is especially valuable for 
solving aerodynamic flow field involving multiple 
scales (i.e. flow field involving regions of very 
low speed in generally high speed flow problems). 
This formulation is also being utilized in 
simulating multi-phase flow such as cavitation in 
liquid flow [36]. 

Numerical Schemes 

It is easy to see that the artificial compressibility 
code resembles that of a compressible flow code. 
For example, to alleviate the difficulty of solving an 
unfactored implicit scheme, indirect methods have 
been devised by many researchers. The alternating 
direction implicit scheme by Beam and Warming 
[37] or Briley and McDonald [38] approximates the 
implicit operator in the unfactored scheme by a 
product of three one-dimensional operators. This 
scheme is unconditionally stable in two dimensions. 
In three dimensions, with the use of numerical 
dissipation, this scheme becomes conditionally 
stable. The ADI scheme introduces factorization 
error that reduces the rate of convergence. However, 
this method introduced perhaps the first implicit 
scheme to CFD, thus viscous flow computation 
became more feasible. With the help of reduced cost 
by the diagonalization of the Jacobian by Pulliam 
and Chaussee [39], this scheme has been used in the 
development of many flow solvers. At NASA 
Ames, the first incompressible Navier-Stokes solver 
in generalized three-dimensional coordinate, INS3D 
code (Kwak et. al [27]), was developed using this 
scheme. The original INS3D code was extensively 
used as a primary design tool for upgraded design of 
the Space Shuttle main engine (SSME) hot gas 
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manifold by Chang et. al [28]. This was probably 
the first major contribution by CFD in rocket 
propulsion system design. 

In artificial compressible formulation, since the 
governing equations are changed into the 
hyperbolic-parabolic type, some of the upwind 
differencing schemes, which have been developed 
for the compressible Euler and Navier-Stokes 
equations [40-43] can be used. These schemes 
were implemented into the flow solver 
successfully using artificial compressibility 
method by Rogers et. al [31]. A variety of 
iterative schemes can be implemented in 
conjunction with the upwind schemes. For 
example, Rogers et. al [31] solved the resulting 
set of numerical equations using a unfactored line 
relaxation scheme similar to that employed by 
MacCormack [44]. 


literature over the years, e.g., problems like square 
cavity, flow over a circular cylinder, duct flow of 
various shapes, backward facing steps, juncture 
flow and tip and wake vortex flow. Even though 
the geometry of these flow problems are relatively 
simple, they offer rich sources of studying flow 
physics and opportunities to prove the validity of 
particular numerical methods and procedures. 
One simple example is shown in Figure 1 where 
the artificial compressibility and pressure 
projection methods are compared to evaluate the 
capability for predicting unsteady phenomena. 
Once the final solutions are obtained, the two 
approaches can produce equivalent results, 
however, the iterative processes for the artificial 
compressibility method and time advancing 
scheme for the pressure projection scheme had to 
be properly handled [48]. 


2.4 Other Methods 

Even though the evolution of incompressible flow 
solver development was explained using examples 
the present authors are directly involved with, a 
large number of research effort has been on-going 
in the incompressible CFD community in general. 
Especially, development along finite element / or 
unstructured grid approach has not been touched 
in this paper, and readers interested in this area 
are referred to a rather recent book by Gresho and 
Sani [9], and a number of articles in the book by 
Hafez [10] (for example, see articles by 
Glowinski [45], Morgan [46], Tezduyar [47] and a 
large number of articles covering numerical 
methods and special topics). Also, a number of 
other methods related to incompressible flow 
computations can be found in the references cited. 



Figure 1. Evolution of stagnation point for an 
impulsively started flat plate: Comparison of artificial 
compressibility (denoted by INS3D-UP), projection 
method (INS3D-FS), finite element solution [49] and 
water tank experiment [50]. 


In general it is the authors opinion that there are 
no major differences any more in basic 
methodologies and tools used among different 
organizations. Differences in robustness, 
computational efficiencies and accuracy usually 
come from how the solution procedures are set up 
in conjunction with geometry and grid handling, 
software engineering and some subtle algorithmic 
aspects not very well presentable in written 
papers. 

3. Computed Examples 

Solutions to some basic problems 

Validation computations of solution algorithms 
and flow solvers have been reported in the 


Validation for real world applications 

Validation of solution procedure for real-world 
applications takes another level of effort because the 
flow field involving complex geometry often 
requires high resolution locally to resolve the 
complicated flow as well as requiring evaluation of 
physical model being implemented. Such validation 
effort is of significant importance to problems like 
submarine with propulsor and turbopump 
subsystems for liquid rocket engine. 

Figure 2 illustrates one such validation using a 
pump impeller. This validation is the first step for 
further application of the simulation procedure to an 
entire turbopump system. Figure 2a shows impeller 
geometry with surface grid where test is conducted 
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for developing an advanced impeller. Figure 2b 
shows some comparison of computed results with 
test data. Even with this simplified geometry, testing 
is not straight forward, which involves variations in 
flow conditions and shroud geometry. Simulation 
involves geometric modeling and grid generation as 
well as approximation related to numerical methods 
and physical modeling. As shown, the computed 
results are not on top of test data unlike the basic 
flow problems cited above. In engineering 
applications, especially, for preliminary design 
applications, one needs to determine what level of 
accuracy is acceptable. More importantly the 
numerical simulation procedure must produce 
consistent results when flow conditions are varied as 
well as for variations in geometry such as different 
angle, length and lean of impeller blades. 



Figure 2a. Geometry of a test impeller. 


Averaged Cm vs. Relative X: R/Rtip = 1.0275 



Figure 2b. Comparison of circumferentially 
averaged meridional velocity at the impeller exit. 

High-end computing for large engineering 
problems 

One of the most significant applications of 
incompressible flow simulation started in the 
early 1980s, when NASA carried out a series of 


upgrades to the Space Shuttle Main Engine 
(SSME) developed in the 1970s. One such effort 
was powerhead redesign. The powerhead, 
considered the backbone of the engine, consists of 
the main injector assembly and pre-burners. 
Partially burned hot gas passes through Hot Gas 
Manifold (HGM) to the main injector assembly. 
During 1983-84, the powerhead redesign was 
undertaken, focusing on two-duct HGM. NASA 
Ames and Rocketdyne collaborated in applying 
the INS3D code for this task. Since the Mach 
number of the partially burned gas was very low, 
the incompressible approximation was a valid 
assumption. The team of researchers successfully 
applied CFD simulation procedure to this task for 
enhancing the performance of the SSME 
powerhead. The new design was developed based 
on two-duct hot gas manifold. The two-duct 
design replaced the previous three-duct engine, 
smoothing the fuel flow, reducing pressure and 
turbulence, and lowering temperatures in the 
engine during operation. This two-duct design 
first flew on the shuttle in July 1995. It 
significantly improved fluid flow in the system, 
thus reducing maintenance and enhancing the 
overall performance of the engine. This 
pioneering work was probably the first major 
application of incompressible CFD to rocket 
propulsion system (see [28] for more detail). 

A large number of incompressible flow software 
has been developed ever since for analyzing 
complex engineering systems. One such example 
is the UNCLE code [36], which has been applied 
to various naval hydrodynamic applications. 
There are, however, many other equally 
significant software in existence. It is not the 
intention of this report to list the inventory of such 
codes. 

A typical process of flow simulation, especially, 
high-fidelity unsteady flow, is illustrated next 
using a turbopump similar to the one used in the 
Shuttle engine (figure 3). 

The goal is to develop a computational framework 
for the design and analysis of an entire liquid 
rocket engine fuel supply system, which is of 
major importance in developing a future rocket 
engine to be used in space launch. A substantial 
computational time reduction for 3D unsteady 
flow simulations is required to reduce the design- 
cycle time of the pumps. Part of this speed up will 
be due to enhancements in computer hardware. 
The remaining portion of the speed-up must be 
contributed by advances in grid generation 
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procedure, flow solution algorithms and by 
efficient parallel implementations. In the present 
example, an overset strategy [51] was adopted to 
generate the grid. To reduce human and 
computing time for generating multiples of pump 
component, scripting [52] has been developed 
thus substantially reducing the grid generation 
portion of the entire simulation time. Figure 4 
illustrate this procedure. 



Figure 3. A typical turbopump for a rocket engine 
showing the geometry and unsteady solution. 



repeated blades, no tip clearance) 
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is based on the explicit message-passing interface 
across MPI groups and is designed for coarse 
grain parallelism. 
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Figure 5. Parallel performance of INS3D-MLP on 
SGI Origin platform. 

The second approach in Multi-Level Parallelism 
(MLP) is obtained by using NAS-MLP routines. 
This approach differs from the MPI/OpenMP 
approach in a fundamental way in that it does not 
use messaging at all. All data communication at 
the coarsest and finest level is accomplished via 
direct memory referencing instructions. This 
approach also provides a simpler mechanism for 
converting legacy code, such as INS3D, than 
MPI. For shared memory MLP, the coarsest level 
parallelism is supplied by spawning off 
independent processes via the standard UNIX 
fork. The boundary data for the overset grid 
system is archived in the shared memory arena by 
each process. Other processes access the data 
from the arena as needed. Figure 5 shows the 
MLP scalability for SSME impeller computations 
using 19.2 Million grid points. 

Biomedical / Biofluid Applications 
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Figure 4. Illustration of automated procedure with 
turbopump scripts for inlet guide vanes. 

Parallel computing strategies vary depending on 
computer architecture such as memory 
arrangement relative to processing units. Two 
approaches have been implemented in the INS3D 
framework: the first approach is hybrid 

MPI/OpenMP and the second one is Multi-Level 
Parallelism (MLP) developed at NASA-Ames 
Research Center. The first approach is obtained 
by using message-passing interface (MPI) for 
inter-zone parallelism, and by using OpenMP 
directives for intra-zone parallelism. INS3D-MPI 


Extension of incompressible flow methods to 
blood flow has been of interest to biomedical 
researchers for many years. For the cardiovascular 
system, the brain and other parts of the human 
body, branching of blood vessels involve 
bifurcations, most of which are non-symmetric. 
Therefore, bifurcation has been a popular subject 
for blood flow simulations. Bifurcation problems 
offer the opportunity to address numerical issues 
involving grid generation as well as to study basic 
fluid dynamic phenomena relevant to blood 
circulation simulation. When the size of blood 
vessel becomes very small as in capillaries, non- 
Newtonian characteristics become significant, 
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thus requiring physical modeling. In addition the 
vessel wall is elastic requiring a structural model 
to account for geometric change depending on 
wall stresses. Early pioneering work on this 
subject was done by Womersley [53]. Later more 
sophisticated model and computations have been 
performed for regular and non-regular bifurcating 
arteries including stenotic vessels (see Berger [54] 
for a review). 

The human circulatory system is like a huge tree 
with many branches of various different sizes. 
Therefore, a lot of computational study has been 
performed using a truncated geometric model. 
One difficulty of simulating a truncated arterial 
system comes from setting proper boundary 
conditions, especially at the downstream 
boundary. To account for large arterial network, 
Quarteroni and his colleagues have a developed 
circuit analog, which has been applied to more 
inclusive circulatory systems modeling such as 
the cardiovascular system and the circle of Willis 
in the brain (see [55] among many publications). 
An alternative way of imposing downstream 
boundary conditions has developed by Tim David 
and his colleagues [56] in their computation of the 
circle of Willis. 

Modeling and simulation of cardiovascular system 
has been the subject of many investigations 
because of its importance in biomedical research. 
Peter Hunter and his colleagues have been 
modeling cardiovascular systems including 
multidisciplinary aspects producing one of the 
most impressive results to date (e.g. see [57]). 
Earlier and perhaps the most elaborate physical 
model of the heart has been pioneered by Peskin 
and his colleagues (see [58] for example). They 
included blood, wall structures and electrical field 
activating heart muscles. 

Another interesting application of blood flow 
simulation is related to artificial devices such as 
artificial heart, ventricular assist device (VAD), 
and heart valves. Since the demand for organs for 
transplantation far exceeds the number of donors, 
the need for artificial devices become increasingly 
high, to be used either as a temporary device or as 
a permanent replacement for natural organ. For 
developing these devices, accurate quantification 
of blood flow plays a crucial role. Thus CFD 
simulation of blood flow in and around these 
artificial devices has become an indispensable 
part of the design. One such example is the recent 
effort of CFD-aided design of DeBakey VAD, 
which enabled human implantation by removing 


thrombus formation and lowering hemolysis to an 
acceptable level for human applications. This 
effort set a new milestone for CFD applications to 
the biomedical area [59], 

All such blood flow computations may be 
regarded as a branch of incompressible CFD. The 
number of CFD applications to blood flow and 
biomedical problems are increasing rapidly, and 
the work cited here represents only small samples 
of the vast amount of work existing and on-going. 

4. Challenges 

Incompressible CFD methodology has been 
advanced along with computational technologies 
in general. Many fluid engineering problems 
could now be simulated, however, mostly at a 
single problem level. For example, it is possible to 
generate solutions to problems like a turbopump 
pump (inducer-impeller-diffuser), a naval vehicle 
at a steady motion including propulsor by a 
model, and a truncated model of the brain or 
heart. To realize the full benefit of CFD, more 
inclusive modeling is desired, such as, systems of 
pump including multiples of pumps and feed line, 
vehicle in maneuver with propulsor, complete or 
more inclusive human circulatory system from 
heart through aorta all the way to brain and kidney 
capillaries. Attempts to solve these have been 
made with some qualitative success. However, the 
predictive capability is still very limited and 
prediction with accurate physics is yet be 
accomplished, which will require inclusion of not 
only fluid dynamic modeling but modeling of 
other quantities like thermal loading, structural 
properties (such as structural behavior of arterial 
wall), and cavitation physics. These computations 
will require not only large computing resources 
but also large data storage and management 
technologies. 

Flow solver codes and software tools have been 
developed to the point that many daily fluid 
engineering problems can now be routinely 
computed. Some of the physical models such as 
for turbulence and transition, however, have not 
been advanced much more than what were 
available in the 70s or 80s. Other models like 
cavitation model are yet to be advanced to 
produce quantitative results for engineering. In 
aerospace design, the most productive aspect of 
CFD applications has been to predict relative 
change among design variations. To push the limit 
of operation and to try bold new ideas, more 
predictive capability will be needed for 
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complicated flows involving transient phenomena, 
separation, tip vortex and cavitation. For example, 
without accurate prediction capability for such 
quantities like cavitation and damaging frequency, 
back flow and rotational stall; CFD cannot be of 
much help to the development of an advance 
turbopump. To make these advances, high-fidelity 
computation using high-end computing facilities 
will still be a “must” despite all the euphoria 
about the PC clusters and grid computing. 

One final note is related to human resources. Even 
though CFD has been advanced remarkably, many 
challenging cases require CFD experts. Computer 
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